# Fig 3. - Map of SBS22a + SBS22b attribution in Romanian and Serbian counties

rm(list = ls())

# Spatial objects
library(geojsonio)
library(sp)
library(sf)
library(raster)
# Data handling
library(dplyr)
library(stringr)
# Plots
library(ggplot2)
library(ggspatial)
library(ggsflabel)
library(ggnewscale)
library(ggpattern)

# Loading individual data
# Residency and duration for all cases along with signature attribution
aa_extended <- readRDS("datasets/aa_extended")

# Estimate of regional effect for SBS22a + SBS22b
sro_regions_estim <- readRDS("datasets/sro_regions_estim") 

# Loading maps
# World map
world_countries <- geojson_read("maps/stanford-kk522dt9425-geojson.json", what = "sp")
# Subset European countries
europe_countries <- world_countries %>% subset(region_un == "Europe")
# Labels for Serbia and Romania only
sro_countries_label <- world_countries %>% subset(name %in% c("Romania","Serbia")) 
kos_label <- world_countries %>% subset(name == "Kosovo")
rm(world_countries)

# Serbian counties
sr_regions <- geojson_read("maps/stanford-hf794fx7640-geojson.json", what = "sp")
# Kosovo 
kos_country <- geojson_read("maps/stanford-pt278dg2451-geojson.json", what = "sp")
# Romanian counties - Grouping Ilfov & Bucuresti
ro_regions <- geojson_read("maps/stanford-rf707cg1788-geojson.json", what = "sp") %>% 
  aggregate(by = "nam")
ro_regions@data <- ro_regions@data %>% mutate(varname_1 = ifelse(nam %in% c("Municipiul Bucuresti", "Ilfov"), "Bucuresti/Ilfov", nam))
ro_regions <- aggregate(ro_regions, by="varname_1")

# Serbian / Romanian country borders: Aggregating all Serbian and Romanian counties
sr_country <- aggregate(sr_regions)
ro_country <- aggregate(ro_regions)
sro_countries <- bind(sr_country, ro_country)

# Kosovo border - Intersection between Kosovo and Serbia
kos <- kos_country %>% st_as_sf() %>% select(geometry) %>% st_cast("LINESTRING")
ser <- sr_country %>% st_as_sf() %>% st_cast("POLYGON")
ser <- ser[1,] %>% st_cast("LINESTRING")
serbia_kosovo_border <- st_intersection(ser,kos) %>% st_cast("POINT") %>% summarise(do_union = FALSE) %>% st_cast("LINESTRING")
rm(kos,ser,sr_country, ro_country)

# Danube and tributaries - Subset from Serbian and Romanian rivers
sr_river <- geojson_read("maps/stanford-wf722gq3793-geojson.json", what = "sp")
ro_river <- geojson_read("maps/stanford-fm546kc1456-geojson.json", what = "sp")
sr_benrivers <- subset(sr_river,nam %in% c("Kolubara","Dunav","Sava","Binacka Morava","Moravica","Velika Morava","Juzna Morava","Zapadna Morava"))
ro_benrivers <- subset(ro_river,str_detect(nam,"Dunarea|Jiu"))
ben_rivers <- bind(sr_benrivers, ro_benrivers)
rm(sr_benrivers, ro_benrivers, sr_river, ro_river)

# Regions of high prevalence for BEN
ben_regions <- geojson_read("maps/map_BEN_regions2.geojson", what = "sp")

# Theme for plot
theme_opts <- list(theme(panel.grid.minor = element_blank(),
                         panel.grid.major = element_blank(),
                         panel.background = element_rect(fill = "lightslategrey"),
                         plot.background = element_rect(fill="white"),
                         panel.border = element_blank(),
                         axis.line = element_blank(),
                         axis.text.x = element_blank(),
                         axis.text.y = element_blank(),
                         axis.ticks = element_blank(),
                         axis.title.x = element_blank(),
                         axis.title.y = element_blank(),
                         plot.title = element_text(size=22)))

# Fig 3
ggplot() +
  # European countries contours & names
  geom_sf(data=st_as_sf(europe_countries), alpha=1, lwd=0.5, fill="ghostwhite")+
  geom_sf_text(data=st_as_sf(europe_countries), aes(label=toupper(name)),alpha=1, size=4.5)+
  # Regional estimates of SBS22a + SBS22b
  geom_sf(data=st_as_sf(sro_regions_estim), aes(fill = AA_sigs_int),  lwd=0)+
  scale_fill_gradient2("Estimate of regional effect on\nSBS22a + SBS22b attribution (Polygons)",low = "#4f81bd", high = "#c0504c", mid = "white", midpoint=0, 
                       guide = guide_colorbar(title.position ="top", barwidth=15, order = 2)) +
  new_scale_fill()+
  # BEN areas
  geom_sf_pattern(data=st_as_sf(ben_regions), alpha=0.9, aes(fill="BEN areas"), pattern='stripe', pattern_colour = "white", pattern_fill = "white", pattern_spacing=.01, pattern_key_scale_factor = 3)+
  scale_fill_manual("",values = c("BEN areas"="sandybrown"), guide = guide_legend(order = 6))+
  new_scale_fill()+
  # Adding contours on top of BEN areas - All countries + Romania and Serbia contours & names
  geom_sf(data=st_as_sf(subset(europe_countries, name %in% c("Bosnia and Herzegovina","Croatia"))), alpha=0, lwd=0.5)+
  geom_sf(data=st_as_sf(sro_countries), alpha=0, lwd=2)+
  geom_sf_label_repel(data=st_as_sf(sro_countries_label),aes(label=toupper(name)), size = 4.5, nudge_x=-2.5,nudge_y=2.2)+
  # Serbia/Kosovo border
  geom_sf(data=serbia_kosovo_border, linetype ="18", lwd=2, size = .1, color="ghostwhite") +
  geom_sf_text(data=st_as_sf(kos_label), aes(label=toupper(name)),alpha=1, size=4.5)+
  # Danube and tributaries
  geom_sf(data=st_as_sf(ben_rivers), aes(colour = "Danube and\ntributaries"),  lwd=1)+
  scale_color_manual("",values = c("Danube and\ntributaries"="steelblue"), guide = guide_legend(order = 7))+
  # Individual locations and values of SBS22a + SBS22b - Zero only
  geom_spatial_point(data = filter(aa_extended, SBS22 == 0 & SBS1536I_cosmic == 0), crs=4326, aes(x = lon_jit, y = lat_jit, size= dur, fill="No attribution of\nSBS22a nor\nSBS22b"),shape=21,alpha=1)+ # removing stroke=2
  scale_fill_manual("",values = c("No attribution of\nSBS22a nor\nSBS22b"="ghostwhite"), guide = guide_legend(order = 4, override.aes = list(size=5)))+
  new_scale_fill()+
  # Individual locations and values of SBS22a + SBS22b - Zero excluded
  geom_spatial_point(data = filter(aa_extended, SBS22 != 0 | SBS1536I_cosmic != 0), crs=4326, aes(x = lon_jit, y = lat_jit, size= dur, fill=AA_sigs),shape=21,alpha=1)+ # removing stroke=2
  scale_fill_viridis_c("SBS22a + SBS22b attribution (log scale)\n(Circles)",direction=-1, trans="log", breaks = c(300,900,2700,8100,24300),
                       guide=guide_colorbar(title.position ="top", barwidth=15,order=5), option="magma")+
  # Scale and north arrow
  annotation_scale(location = "bl", width_hint = 0.2, text_cex = 1) +
  annotation_north_arrow(location = "bl", which_north = "true",
                         pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
                         style = north_arrow_fancy_orienteering) +
  coord_sf(crs=4326,xlim = c(14.1, 29.8), ylim = c(42.2,48.3), expand = FALSE) +
  theme(legend.position = "bottom", legend.key = element_rect(fill="white"),
        legend.key.height = unit(0.5, "cm"), legend.text = element_text(size=12), legend.title = element_text(size=12)) +
  guides(size=guide_legend(title.position="top", order=1,nrow=2,byrow=T)) +
  labs(size="Duration of residency\n(in years)")+
  theme_opts
# Output PDF in 9x14
